# libraries
library(plyr)
library(tidyverse)
## Warning: Installed Rcpp (0.12.12) different from Rcpp used to build dplyr (0.12.11).
## Please reinstall dplyr to avoid random crashes or undefined behavior.
## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr
## Conflicts with tidy packages ----------------------------------------------
## arrange():   dplyr, plyr
## compact():   purrr, plyr
## count():     dplyr, plyr
## failwith():  dplyr, plyr
## filter():    dplyr, stats
## id():        dplyr, plyr
## lag():       dplyr, stats
## mutate():    dplyr, plyr
## rename():    dplyr, plyr
## summarise(): dplyr, plyr
## summarize(): dplyr, plyr
library(ggplot2)
# Read files from one folder and bind them together into one df

data_directory_path = "testData"

filenames <- list.files(data_directory_path, pattern="*.csv", full.names=TRUE)
list_of_dataframes <- lapply(filenames, read.csv)
busData <- rbind.fill(list_of_dataframes)

Primeiramente, vamos ver as variáveis que existem nos dados:

str(busData)
## 'data.frame':    190858 obs. of  37 variables:
##  $ route                    : int  203 203 203 203 203 203 203 203 203 203 ...
##  $ tripNumOrig              : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ shapeId                  : int  3855 3855 3855 3855 3855 3855 3855 3855 3855 3855 ...
##  $ shapeSequence            : int  4649120 4649137 4649154 4649184 4649211 4649239 4649258 4649278 4649300 4649329 ...
##  $ shapeLatOrig             : num  -25.5 -25.5 -25.5 -25.5 -25.5 ...
##  $ shapeLonOrig             : num  -49.3 -49.3 -49.3 -49.3 -49.3 ...
##  $ distanceTraveledShapeOrig: num  485 847 1256 1778 2243 ...
##  $ busCode                  : Factor w/ 544 levels "","BA013","BA124",..: 46 46 46 46 46 46 46 46 46 46 ...
##  $ gpsPointId               : logi  NA NA NA NA NA NA ...
##  $ gpsLat                   : num  -25.5 -25.5 -25.5 -25.5 -25.5 ...
##  $ gpsLon                   : num  -49.3 -49.3 -49.3 -49.3 -49.3 ...
##  $ distanceToShapePoint     : num  7.75 24.43 2.65 5.63 7.77 ...
##  $ timestampOrig            : Factor w/ 60313 levels "","05:11:15",..: 28 30 36 45 52 60 68 82 89 100 ...
##  $ busStopIdOrig            : int  25739 25737 25540 26197 25924 25925 26009 25523 26010 25733 ...
##  $ problem                  : Factor w/ 4 levels "BETWEEN","NO_PROBLEM",..: 2 2 2 2 2 2 1 2 2 2 ...
##  $ date                     : Factor w/ 7 levels "2017-02-01","2017-02-02",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ busStopIdDest            : int  25737 25540 26197 25924 25925 26009 25523 26010 25733 25732 ...
##  $ timestampDest            : Factor w/ 60335 levels "","05:12:50",..: 27 33 40 47 55 63 76 83 94 106 ...
##  $ tripNumDest              : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ shapeLatDest             : num  -25.5 -25.5 -25.5 -25.5 -25.5 ...
##  $ shapeLonDest             : num  -49.3 -49.3 -49.3 -49.3 -49.3 ...
##  $ distanceTraveledShapeDest: num  847 1256 1778 2243 2821 ...
##  $ duration                 : int  41 101 100 85 72 55 94 58 70 92 ...
##  $ distance                 : num  363 409 522 466 578 ...
##  $ hourOrig                 : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ hourDest                 : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ isRushOrig               : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ isRushDest               : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ periodOrig               : Factor w/ 3 levels "afternoon","morning",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ periodDest               : Factor w/ 3 levels "afternoon","morning",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ weekDay                  : Factor w/ 7 levels "Fri","Mon","Sat",..: 7 7 7 7 7 7 7 7 7 7 ...
##  $ weekOfYear               : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ dayOfMonth               : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ month                    : int  2 2 2 2 2 2 2 2 2 2 ...
##  $ isHoliday                : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ isWeekend                : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ isRegularDay             : int  1 1 1 1 1 1 1 1 1 1 ...

Podemos ver que existem 37 variáveis, onde temos: * route : rota do ônibus * tripNumOrig : o número da viagem do ônibus na origem * shapeId : o id do shape do ônibus * shapeSequence : o ponto do shape que a viagem está atualmente * shapeLatOrig : a latitude da origem do shape * shapeLonOrig : a longitude da origem do shape * distanceTraveledShapeOrig: a distancia que o ônibus já viajou na origem. * busCode : o código do ônibus * gpsPointId : o id do ponto do gps * gpsLat : a latitude do ponto do gps * gpsLon : a longitude do ponto do gps * distanceToShapePoint : a distância do gps para o ponto do shape casado * timestampOrig : o horário no ponto de origem * busStopIdOrig : a parada do ônibus no ponto de origem * problem : o tipo de problema encontrado pelo BULMA que pode ser NO_PROBLEM, quando não houver problema, NO_SHAPE, quando nenhuma trip foi formada (gps sem shape na base de dados ou com pontos insuficientes para formar trip), TRIP_PROBLEM, quando não foi possível identificar ponto final/ponto inicial, OUTLIER_POINT, quando os pontos excederam o threshold de distância para o shape. * date : a data da viagem * busStopIdDest : o ponto de ônibus do destino * timestampDest : o horário no ponto de destino * tripNumDest : o número da viagem do ônibus no destino * shapeLatDest : a latitude do ponto do shape destino * shapeLonDest : a longitude do ponto do shape destino * distanceTraveledShapeDest: a distancia que o ônibus já viajou no destino * duration : a duração entre a parada origem e destino * distance : a distancia entre a parada origem e destino * hourOrig : a hora de origem * hourDest : a hora de destino * isRushOrig : indicativo se está em horário de pico na origem * isRushDest : indicativo se está me horário de pico no destino * periodOrig : o período na origem (manhã, tarde, noite) * periodDest : o período no destino (manhã, tarde, noite) * weekDay : o dia da semana * weekOfYear : a semana do ano * dayOfMonth : o dia do mês * month : o mês * isHoliday : indicativo se está em mês de férias * isWeekend : indicativo se é final de semana * isRegularDay : indicativo se está entre a terça e quinta inclusive

Análises posssíveis:

  1. um panorama das distâncias entre as paradas
  2. um panorama do número de shapes por rota
  3. um panorama do número de paradas por shape
  4. um panorama do tamanho dos shapes
  5. Relação entre duração entre uma parada e outra e a distância entre essas paradas
  6. panorama do periodo do dia e duração da viagem entre duas paradas
  7. Panorama entre o dia da semana e a duração da viagem entre duas paradas

analises que queremos:

  1. Relação entre a hora do dia e duração de uma viagem entre duas paradas
  2. Relação entre o dia da semana e a duração de uma viagem entre duas paradas
  3. Relação entre a área da cidade e a duração de uma viagem entre duas paradas

Hora do dia

Agora vamos plotar um gráfico de duração de viagens entre paradas por horas do dia.

p <- ggplot(busData, aes(x = as.factor(hourOrig), y = duration, colour = problem))
p + geom_jitter(aes(alpha = 0.5))
## Warning: Removed 5745 rows containing missing values (geom_point).

Problemas encontrados:

  1. Algumas durações estão negativas (pouquíssimos casos 8 em 190 mil da amostra).
  2. Algumas durações estão muito longas. Muitos desses casos ocorrem quando o número de viagem da origem é diferente do número de viagem do destino.

Há um problema de cálculo no preprocessamento quando a parada de origem é a última parada do shape e a parada de destino é a primeira parada do shape. Filtrando esses casos obtemos o gŕafico abaixo.

temp <- busData %>% filter(distanceTraveledShapeDest > 0)

p <- ggplot(temp, aes(x = as.factor(hourOrig), y = duration, colour = problem))
p + geom_jitter(aes(alpha = 0.5))
## Warning: Removed 3401 rows containing missing values (geom_point).

Há casos em que o número da viagem do ponto de origem e o número da viagem do ponto de destino não são iguais (nem consecutivos?). Filtrando esses casos obtemos o gŕafico abaixo.

temp2 <- temp %>% filter(abs(tripNumOrig - tripNumDest) == 0)

p <- ggplot(temp2, aes(x = as.factor(hourOrig), y = duration, colour = problem))
p + geom_jitter(aes(alpha = 0.5))
## Warning: Removed 1757 rows containing missing values (geom_point).

Não é interessante que haja casos em que a parada de origem é igual à parada de destino. Vamos filtrar esses casos e visualizar como ficam os dados.

temp3 <- temp2 %>% filter(busStopIdOrig != busStopIdDest)

p <- ggplot(temp3, aes(x = as.factor(hourOrig), y = duration, colour = problem))
p + geom_jitter(aes(alpha = 0.5))
## Warning: Removed 614 rows containing missing values (geom_point).

Não mudou muita coisa. Há casos em que simplesmente o ônibus demorou muito para ir de uma parada a outra e não se sabe exatamente o porquê.

Podemos calcular a velocidade média do ônibus para ir de uma parada a outra. O resultado está no gráfico abaixo.

temp3 <- temp3 %>% mutate(velocityKmh = ((distance / duration) * 3.6))

p <- ggplot(temp3, aes(x = "velocidade", y = velocityKmh, colour = problem))
p + geom_jitter(aes(alpha = 0.5))
## Warning: Removed 614 rows containing missing values (geom_point).

Às vezes a velocidade é infinita (quando a duração é 0). Há muitas velocidades que não fazem sentido para um ônibus. A velocidade máxima para ônibus em trechos urbanos é de 90 km/h.

# Filtrando velocidades irreais (muito altas)
temp4 <- temp3 %>% filter(velocityKmh <= 90)

# filtrando viagens sem problemas e interpoladas
temp5 <- temp4 %>% filter(problem %in% c("NO_PROBLEM", "BETWEEN"))

p <- ggplot(temp5, aes(x = "velocidade", y = velocityKmh))
p + geom_jitter(aes(alpha = 0.5, colour = problem))

p + geom_boxplot(aes(colour="red"))

É possível perceber que cerca de 75% das observações tem velocidade abaixo de 35 km/h, sendo a mediana próxima de 25 km/k.

p <- ggplot(temp5, aes(x = as.factor(hourOrig), y = duration, colour = problem))
p + geom_jitter(aes(alpha = 0.5))

É notório que existem alguns outliers na visualização acima.

Vamos ver um ponto de corte aceitável para a duração:

percentiles <- data.frame(perc = seq(0, 1, 0.0001))

percentiles <- percentiles %>% 
  mutate(
    value = quantile(temp5$duration, perc),
    diff = value - lag(value)     
         )

ggplot(percentiles, aes(x = perc, y = value)) +
  geom_line()

Podemos ver que os dados estão bem concentrados até certo ponto, então vejamos as últimas observações para descobrirmos um bom ponto de corte de outliers:

percentiles %>% 
  arrange(-perc) %>% 
  head(50)
##      perc      value       diff
## 1  1.0000 52453.0000 45782.3379
## 2  0.9999  6670.6621  2362.6093
## 3  0.9998  4308.0528  1134.1679
## 4  0.9997  3173.8849   510.7117
## 5  0.9996  2663.1732   242.4942
## 6  0.9995  2420.6790   191.6790
## 7  0.9994  2229.0000   319.0000
## 8  0.9993  1910.0000   147.8944
## 9  0.9992  1762.1056   144.6304
## 10 0.9991  1617.4752   133.4942
## 11 0.9990  1483.9810    93.1584
## 12 0.9989  1390.8226    82.7434
## 13 0.9988  1308.0792    67.9113
## 14 0.9987  1240.1679    63.1679
## 15 0.9986  1177.0000    41.0000
## 16 0.9985  1136.0000    47.0000
## 17 0.9984  1089.0000    31.9567
## 18 0.9983  1057.0433    37.4245
## 19 0.9982  1019.6188    28.6811
## 20 0.9981   990.9377    14.9377
## 21 0.9980   976.0000    21.6114
## 22 0.9979   954.3886    19.3886
## 23 0.9978   935.0000    13.5491
## 24 0.9977   921.4509    22.3717
## 25 0.9976   899.0792    17.0792
## 26 0.9975   882.0000    10.6642
## 27 0.9974   871.3358    11.3358
## 28 0.9973   860.0000    18.2228
## 29 0.9972   841.7772    23.7772
## 30 0.9971   818.0000    13.0000
## 31 0.9970   805.0000     8.5227
## 32 0.9969   796.4773    22.3717
## 33 0.9968   774.1056    10.1056
## 34 0.9967   764.0000    12.2756
## 35 0.9966   751.7244    13.7339
## 36 0.9965   737.9905     6.7529
## 37 0.9964   731.2376     9.7434
## 38 0.9963   721.4942     8.4942
## 39 0.9962   713.0000     6.0000
## 40 0.9961   707.0000     5.0000
## 41 0.9960   702.0000     4.0000
## 42 0.9959   698.0000     6.0000
## 43 0.9958   692.0000     6.9831
## 44 0.9957   685.0169     8.0169
## 45 0.9956   677.0000     6.0000
## 46 0.9955   671.0000     5.0000
## 47 0.9954   666.0000     5.0000
## 48 0.9953   661.0000     4.8416
## 49 0.9952   656.1584     4.1584
## 50 0.9951   652.0000     3.0000

Podemos ver que a diferença entre os valores dos percentis, os quais aumentam em 0.001 a cada linha, começam a aumentar repentinamente a partir do percentil 0.9984. Com isso, vamos estabelecer que a duração é um outlier quando for maior que 1100 segundos.

Vamos verificar se este ponto de corte de duração faz sentido em relação ao geral:

threshold <- 1100
temp5 %>% 
  ggplot(aes(y = duration, x = problem, colour=problem)) +
  geom_jitter(aes(alpha = 0.5)) +
  geom_violin() +
  geom_hline(yintercept = threshold)

Podemos perceber que no geral, as durações acima do ponto de corte estão bem dispersas, enquanto as que estão abaixo são mais concentradas.

Com isso, vamos filtrar as observações que tem a duração menor que o limiar.

temp6 <- temp5 %>% filter(duration <= threshold)

p <- ggplot(temp6, aes(x = as.factor(hourOrig), y = duration))
p + geom_jitter(aes(alpha = 0.5, colour = problem))

p + geom_boxplot()

Vemos que alguns pontos ainda estão dispersos, mas estão bem mais concentrados do que antes do corte.

Além disso, podemos ver que entre as 7 e 18 horas a duração das viagens entre paradas são bastante parecidas, enquanto as demais tem valores mais baixos.

write.csv(x = temp6, file = "dadostratados.csv", row.names = F)